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Abstract.  We  solve  statistical  moment  differential  equations  (MDEs)  for  immiscible  flow  in 
porous  media  in  the  limit  of  zero  capillary  pressure,  with  application  to  secondary  oil  recovery. 
Closure  is  achieved  by  Taylor  expansion  of  the  fractional  flow  function  and  a  perturbation  argument. 
Previous  results  in  1-D  are  extended  to  2-D,  in  which  a  bimodal  profile  is  less  evident. 

Mean  and  variance  of  (water)  saturation  exhibit  a  bimodal  character;  two  shocks  replace  the 
single  shock  front  evident  in  the  classical  Buckley-Leverett  saturation  profile.  Comparison  to  Monte 
Carlo  simulations  (MCS)  shows  that  the  MDE  approach  gives  a  good  approximation  to  total  oil 
production.  For  such  integrated  or  averaged  quantities,  or  where  a  rough  approximation  of  the 
location  and  magnitude  of  uncertainty  is  sufficient,  MDEs  may  be  substantially  more  efficient  than 
MCS. 
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1.  Introduction.  Stochastic  representations  of  subsurface  geologic  properties 
have  become  commonplace  due  to  the  difficulty  in  complete  and  certain  characteriza¬ 
tion  of  these  properties.  This  leads  to  uncertainty  in  flow  profiles  in  such  porous  media, 
so  that  statistical  description  of  outcomes  is  appropriate.  Additionally,  from  the  per¬ 
spective  of  practical  macroscopic  field-scale  models,  microscopic  heterogeneities  and 
flows  can  be  viewed  as  random  processes  to  be  upscaled.  In  two-phase  flow  on  a  field 
scale,  we  are  primarily  interested  in  mean  behavior  and  a  measure  of  uncertainty  in 
the  saturation  of  each  phase,  which  can  be  a  measure  of  “macrodispersion”  in  certain 
upscaling  contexts. 

A  “zeroth-order”  model  of  mean  flow  with  averages  of  geologic  properties  ignores 
correlations  between  flow  variables.  Monte  Carlo  simulations  (MCS)  of  many  real¬ 
izations  of  geologic  properties  to  estimate  moments  require  much  computation  time 
and  careful  sampling  techniques  [10,11,29].  Macrodispersion  theory  in  contaminant 
transport  captures  an  approximate  effect  of  fluctuations  by  modeling  covariance  func¬ 
tions  in  an  equation  for  the  mean  concentration  [5, 12].  This  theory  has  a  long  history 
in  subsurface  contaminant  transport  and  is  closely  related  to  eddy  diffusion  models 
of  turbulence  [15,  p.  358  ff] . 

Rather  than  approximate  covariance  functions  as  macrodispersive  terms,  which 
requires  neglecting  fluctuations  in  second  moments,  we  solve  PDEs  for  the  covari¬ 
ance  functions  and  the  mean  of  saturation,  extending  previous  results  to  2-D.  These 
moment  differential  equations  (MDEs)  are  derived  using  a  modified  perturbation  ex¬ 
pansion  as  described  in  [17].  The  MDEs  allow  direct  approximation  of  the  local  mean 
and  covariance  functions  for  general  boundary  conditions  and  general,  nonstationary 
stochastic  geology  [33]. 
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1.1.  Applications.  A  statistical  description  of  subsurface  flow  is  of  particular 
interest  for  secondary  oil  recovery.  In  the  advection  equation,  the  principal  difficulty  is 
a  non-convex  nonlinear  flux  function  that  leads  to  discontinuous  solutions.  Standard 
pressure-saturation  equations  for  2-D  horizontal  flow  of  two  incompressible  immiscible 
fluids  in  porous  media  in  the  limit  of  vanishing  capillary  pressure  ( pc  =  0)  are 

</>v(x)  =  —  K(x)  V/i(x),  V  •  v(x)  =  0,  (1.1) 

<9ts(x,  t)  +  V  •  [/(s(x,  i))v(x)]  =  0.  (1.2) 

These  are  considered  valid  from  laboratory  (centimeters)  to  field  scales  of  reservoir 
depth  (10-100  meters)  and  length  (100-10,000  meters).  Hydraulic  conductivity  K 
may  be  an  anisotropic  tensor;  here,  for  simplicity,  it  will  be  an  isotropic  scalar  K. 
Assume  also  that  K  depends  weakly  on  (water)  saturation  s  [1,3].  Apply  (1.1)-(1.2) 
to  the  flow  of  oil  and  water,  for  arbitrary  fluid  mobilities.  Denote  the  total  velocity, 
a  scaled  total  volumetric  flux  of  both  fluids,  by  v,  hydraulic  water  head  by  h ,  and 
porosity  (assumed  constant)  by  f>.  The  fractional  flow  function  f(s),  for  pc  =  0, 
represents  the  fraction  of  v  due  to  water.  It  is  typically  S-shaped  as  shown  in  Fig. 
1.1;  in  §4.1  we  use  a  form  arising  from  quadratic  relative  permeabilities  (see  [1]). 
However,  our  method  does  not  depend  on  any  such  specific  choice  of  /(s). 


Fig.  1.1.  The  nonconvex  fractional  flow  function. 

Capillary  pressure  regularizes  sharp  fronts  caused  by  the  nonlinear  advection 
term.  To  obtain  a  linear  approximation  to  this  effect,  add  e£>V2s(x,<)  with  cd  >  0 
to  the  right  side  of  (1.2).  Letting  >  0  defines  the  vanishing-viscosity  solution  [27], 
which  is  the  one  we  seek. 

As  is  standard  in  subsurface  applications,  let  Y  =  In  K  be  a  random  field  with 
prescribed  mean  and  covariance  functions;  e.g.,  it  is  often  claimed  that  Y  is  multivari¬ 
ate  Gaussian,  based  on  empirical  observations  [8]  (again  our  method  does  not  depend 
on  such  a  limiting  assumption).  Through  (1.1) — (1.2),  v  and  s  are  thus  random  fields. 
No  other  underlying  sources  of  uncertainty  are  considered  in  this  study. 

Under  the  assumptions  stated  above,  with  steady  boundary  conditions,  and  ne¬ 
glecting  the  weak  dependence  of  K  on  s,  a  steady  v  can  be  determined  from  (1.1). 
We  evolve  s  from  the  stochastic  PDE  (1.2),  assuming  known  statistics  of  v.  Moments 
of  v  and  h  can  be  estimated  from  established  theory  ([31-33,36]  use  MDEs).  We 
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combine  analytical  and  numerical  techniques  to  model  the  propagation  of  uncertainty 
from  an  underlying  random  field  V(x),  through  v(x)  to  the  solution  s(x,t). 

1.2.  Previous  work.  Existing  work  on  MDEs  focuses  mostly  on  advection  equa¬ 
tions  with  linear  flux  functions  [13]  and  some  nonlinear  subsurface  flow  equations  of 
a  form  different  from  (1.2)  [2,  28,  31,  32,  36]. 

Langlo  and  Espedal  [19, 20]  present  a  macrodispersion  approach  for  the  stochastic 
version  of  (1-2).  The  flux  function  is  expanded  in  a  Taylor  series,  and  high-order  terms 
are  neglected;  then  standard  techniques  represent  macrodispersivity  as  a  function  of 
flow  velocity  covariance.  Zhang,  Tchelepi,  and  Li  take  advantage  of  the  steady  velocity 
field,  and  transform  2-D  flow  to  1-D  Lagrangian  flow  along  streamlines  [34,  35] .  They 
formulate  integral  equations  for  moments  from  ensemble  averages  over  the  streamlines 
rather  than  a  system  of  MDEs.  Variations  on  the  streamline  approach  are  currently 
popular  in  stochastic  subsurface  hydrology  [4,  6,  7,  25, 26]. 

An  Eulerian  MDE  approach  has  been  successful  for  single-  and  multiphase  pres¬ 
sure  and  velocity  equations  [31],  and  a  natural  next  step  is  to  extend  the  theory  from 
flux  equations  to  transport  equations  (as  in  [13])  and  (1.2).  The  Eulerian  framework 
differs  from  streamlines  not  only  in  formulation  but  also  in  that  the  MDEs  need  no 
velocity-distribution  assumption  and  extension  to  transient  velocity  fields  is  relatively 
straightforward.  The  appeal  of  MDEs  relative  to  macrodispersion  approaches  is  that 
covariances  are  computed  directly,  so  that  it  is  not  necessary  to  approximate  them 
by  solving  an  approximate  fluctuation  equation  that  neglects  second-moment  fluctu¬ 
ations.  In  [17]  and  [18],  we  derived  and  solved  second-order  MDEs  for  (1.2)  in  1-D, 
using  three  different  perturbation  approaches.  One  of  these  approaches,  a  modified 
expansion  described  in  [17],  is  applied  here  in  2-D. 

Equations  are  presented  in  §2.  In  §3,  we  discuss  classification  of  MDEs  and  the 
fundamental  difference  between  2-D  and  1-D  MDEs.  Numerical  solutions  are  given 
in  §4  and  are  compared  to  Monte  Carlo  simulations. 

2.  Moment  Equations.  We  use  a  modified  perturbation  expansion,  described 
in  detail  in  [17],  to  derive  statistical  MDEs.  Closure  implicitly  depends  upon  the 
assumption  that  random  fluctuations  in  Y  =  In  K  are  small  (cry  <C  1).  Moments  of  h 
and  v  are  assumed  known  from  (1.1)  and  moments  of  Y. 

The  saturation  equation  (1.2),  with  initial  data  s(x,0)  =  c/(x),  is 


dts  +  dXi{f{s)v.i)  =  0  (2.1) 

(Einstein  summation  convention  assumed) .  We  assume  that  g  is  known  with  certainty. 
Solutions  are  defined  in  terms  of  vanishing  viscosity  as  in  §1.1;  henceforth,  this  is 
tacitly  understood. 

For  examples  of  commonly  used  methods  for  deriving  moment  equations  with 
applications  to  subsurface  flow  and  transport,  see  [5,9,12,13,23,31].  These  meth¬ 
ods  originated  in  the  correlation  equations  of  turbulence  models.  Here  we  apply  a 
modification  to  a  standard  asymptotic  expansion. 

Let  (•)  denote  the  expectation  operator,  defined  by 

(ip)  =  j  ip(u)  dP(uj)  (2.2) 

Jn 

for  any  integrable  function  ip  :  — >  R.  on  the  sample  space  O  with  probability  measure 

P.  We  omit  reference  to  to  in  what  follows. 
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The  random  field  Y  is  decomposed  into  deterministic  mean  plus  random  fluctu¬ 
ation:  Y  =  (Y)  +  5Y.  Each  field  dependent  on  Y  is  represented  by  a  formal  power 
expansion  in  an  as-yet  unknown  parameter  e: 

OO  OO  OO 

h  =  e"Mx)>  v  =  ^2  e”v„(x),  s  =  ^  ensn(x,t).  (2.3) 

n=0  n= 0  n= 0 

The  expansion  parameter  e  =  ay  is  shown  to  be  appropriate  within  the  context  of 
the  velocity  and  head  equations  (1.1)  [5,  pp.  184-190],  [33].  For  example,  for  single¬ 
phase,  stationary  uniform  mean  flow  on  an  infinite  domain  in  1-D,  vq  is  a  deterministic 
scalar  and  cr2  is  shown  to  be  approximated  by  e2  (r>2)  =  Vq a'y .  In  2-D,  we  denote  the 
components  of  vn  by  (i’i)n  for  *  =  1,2. 

In  the  following,  the  term  “order”  applies  to  the  power  of  e  rather  than  the  order 
of  the  statistical  moment.  In  special  cases,  we  may  use  the  term  “order”  in  the 
latter  sense  and  clearly  differentiate  from  order  in  e.  Thus,  a  “second  moment”  is  a 
covariance,  but  “second-order  moment”  is  any  moment  that  is  approximated  to  order 
e2.  Recall  that  we  only  need  the  decompositions  of  v  and  s  here.  The  fractional  flow 
function  is  expanded  in  a  Taylor  series  around  s  =  So  +  e  (si)  +  e2  ($2)  (a  second-order 
approximation  to  the  mean  saturation;  we  see  below  that  (so)  =  so): 

/(s)  =  /(®)  +  f(s)(s  ~  s)  +  - f"(s)(s  -  s)2  d  ,  (2.4) 

where 

OO 

s  —  s  =  Ssq  +  eSsi  +  e2Ss2  +  e™s„(x,  t). 

n— 3 

The  expansion  center  s  is  our  approximation  to  the  mean  saturation,  so  we  retain  s 
as  the  argument  of  /  in  the  following  modification  to  the  standard  equations  of  order 
0, 1,  and  2.  The  modified  equations  for  sq  and  Si,  obtained  by  substituting  (2.3)  and 
(2.4)  into  (2.1)  and  collecting  the  powers  e°  and  e,  are  given  by 

dts0  +  dXi[f(s)(vi)0]  =  0,  (2.5) 

and 

dtS!  +  dXi  f(s)(v.i)i  + f(s)5si(vi)0  =0.  (2.6) 

Note  that  (so)  =  so,  because  all  other  terms  in  the  equation  for  so  are  deterministic. 
We  show  that  (si)  =  0.  By  definition,  (<5si)  =  0.  Also,  (vi)  =  0  by  a  standard 
result  from  the  theory  for  moment  equations  of  multidimensional  single-phase  flow 
[36].  Applying  the  expectation  operator  (•)  to  the  first-order  equation  (2.6)  gives 

0  =  dt  (si)  +  dXi  f(s)  ((wj)i)  +  f(s)  (Ssx)  (vi)0  =dt(s1).  (2.7) 

The  initial  data  are  deterministic;  therefore,  (si)  is  initially  zero,  and  (2.7)  implies 
that  (si)  remains  zero  for  all  time.  The  mean  approximation  is  then  s  =  So  +  e2  (S2). 
The  correction  term  S2  satisfies 

dtS2  +  dXi  f(s)(vi) 2  +  f'(s)6si(vi)i  +  f'(s)5s2(vi) 0  +  ^/"(s)i5si2(^)o  =  0. 
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Apply  the  expectation  operator  and  note  that  ( <5si 2)  =  (siSi)  —  (si)  (si)  =  (s i2), 
(<5s2)  =  0,  and  Vi  =  <5vi  implies  ($si(t>j)i)  =  (<5si(^)i)  =  (si(uj)i)  -  (si)  ((v»)  1)  = 
to  obtain  the  following  correction  to  (2.5): 

dt(s2)  +  dXi  f(s)  ((vi)2)  A  f'(s)  (si(ui)i)  +  ^/"(s)  (si2)  (vi)0  =0.  (2.8) 

Define  v  =  Vo  +  e2  (v2),  which  is  consistent  with  the  second-order  theory  for  single¬ 
phase  flow.  Adding  the  zeroth-order  equation  (2.5)  to  e2x(2.8)  gives  an  equation  for 
the  mean  saturation: 

dts  +  dXi  f(s)vi  +  e2f(s)  (si(ui)i)  +  ^-/"(s)  (si2)  (wj)o  =0.  (2.9) 

We  derive  equations  for  the  unknown  covariance  functions  (si(uj)i)  and  (s2) 
using  the  following  additional  notation.  The  independent  variables  are  x  and  t  except 
where  noted,  and  -|y  denotes  the  replacement  of  x  by  some  y  different  from  x.  It  is 
convenient  and  useful  to  derive  equations  for  the  more  general  two-point  covariances 
(si(uj)i|y)  and  (sisi|y)  rather  than  for  one-point  covariances. 

To  obtain  an  equation  for  (si(uj)i|y),  multiply  (2.6)  by  (uj)i(y),  and  apply  (•). 
This  results  in 

dt  (si(ui)i|y)  +  dXi  f(s)  ((u»)i(ui)i|y)  +  f'(s)(vi)o  (si(uj)i|y)  =0.  (2.10) 

Similarly,  multiply  (2.6)  by  si(y,i)  and  use  the  identity1  Si|y9tSi  =  3t(siSi|y)  — 
si^tSi|y,  then  use  (2.6)  with  y  in  place  of  x,  yielding  this  equation  for  the  two-point 
saturation  covariance: 

(^]Al|y)  A  &Xi  /(^)  (^lly(^i)l)  A  f  (s)(Ui)Q  (Sl-Si|y) 

/(S|y)  (si(uj)i|y)  +  /  (s|y)(uj)o|y  (si|ySi)  =0.  (2-11) 

Define  ca„4(x,y,t)  =  e2  (si(uj)i|y) ,  cs  =  e2(siSi|y),  cViV.  =  e2  ((w*)i(uj)i ly)> 
(s)  |y  =  (s)  (y,i),  (7 SVi  (x,  t)  =  cSVi(x,xi,t),  and  <r2(x,f)  =  cs(x,x,t).  Let  a  caret  over 
a  variable  denote  the  mapping  ?/>(x,  y,  t)  =  ip(y,  x,  t)  for  any  function  ip.  The  resulting 
system  is 

dts  A  dXi  f(s)vi  +  f(s)aaVi  +  ^/,,(S)(ui)0  <J2S  =  0,  (2.12a) 

dtcSVj+dXi  f(s)cViV.+ f{s)(vi)0cSVj  =0,  (2.12b) 

dtcs  +  dXi  f(s)cSVi  +  f'{s)(vi)0cs  +  dVi  f(s)cSVi  + f(s)(vi)Q\yCs  =0;  (2.12c) 

i,j  =  1,2. 

Recall  that  velocity  moments  (vj)o,  u*,  and  cViV.  are  assumed  known.  The  flux  in 
(2.12a)  consists  of  a  nonlinear  advective  mean  flux  term  and  two  covariance  terms. 
In  turbulence  applications,  terms  such  as  cSVj  often  are  referred  to  as  transport  by 
fluctuations  [22].  Both  (2.12b)  and  (2.12c)  have  advective  flux  terms,  are  coupled  to 

1  The  identity  is  not  valid  in  a  strong  sense  for  discontinuous  solutions.  Recall,  however,  that  we 
define  solutions  in  terms  of  the  (smooth)  viscous  solution,  in  the  limit  ep  — >  0. 
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the  mean  equation  (2.12a),  and  are  first-order  in  e2  =  ay.  This  is  consistent  with  the 
approximation  s,  which  is  second-order  in  ay.  However,  due  to  the  choice  of  s,  higher- 
order  effects  are  included  in  (2.12).  We  demonstrate  this  fact  for  the  case  f(s)  =  s, 
where  the  mean  advective  flux  term  in  (2.12a)  is  dXi  [su,:],  or  dXi [(so  +  e2  (s2))((fi)o  + 
e2  ((ui)2))].  The  term  e4dXi  [(S2)  ((v*)2)]  in  this  expression  is  a  partial  fourth-order 
correction. 

We  refer  to  this  problem  as  2-D  although  the  covariance  functions  involve  four 
spatial  variables  and  the  saturation  covariance  has  fluxes  in  all  directions.  The  four 
variables  represent  two  different  points  in  the  same  2-D  domain. 

To  formally  define  the  vanishing-viscosity  solution,  first  add  the  deterministic 
term  eodx.s  to  (2.1)  and  carry  out  the  expansion  and  derivation  above.  This  adds 
diffusion  terms  end^.s,  eodx.cSVj,  and  e_D<92.cs  to  the  right-hand  sides  of  (2.12a), 
(2.12b),  and  (2.12c),  respectively.  Then  find  solutions  in  the  limit  eo  — »■  0.  In  practice 
we  solve  the  system  (2.12)  directly. 

3.  Classification.  Classification  is  of  interest  for  several  reasons.  Macrodisper¬ 
sion  theory  produces  a  parabolic  equation  for  the  mean,  with  a  macrodispersivity  that 
depends  on  cViVj  and  (s).  In  contrast,  we  have  shown  that  1-D  MDEs  are  hyperbolic 
[18].  To  achieve  this  result,  we  exploited  special  structure  to  reduce  the  equations. 
Namely,  1-D  velocity  has  an  infinite  correlation  length  because  it  is  a  spatially  constant 
random  variable.  This  led  to  the  covariance  relationship  cs  cv  =  csv  csv.  Velocity 
correlations  have  a  finite  length  scale  in  2-D.  We  argue  that  this  fact  prohibits  classifi¬ 
cation  and  suggests  a  transition  to  parabolic  behavior  as  described  by  macrodispersion 
theory. 

We  expect  the  equations  to  be  nearly  hyperbolic  since  the  original  deterministic 
PDE  (1.2)  is  hyperbolic.  However,  <jSVi  and  <r2,  rather  than  cSVi  and  cs,  appear 
in  the  mean  equation.  This  mix  of  one-point  and  two-point  covariance  functions 
prevents  us  from  immediately  treating  the  MDEs  as  classically  hyperbolic.  Also, 
independent  variables  x  and  y  are  permuted  in  s  and  cSVj  in  (2.12c).  In  general, 
s | y  =  s(y,  t)  yf  5(x,t),  and  csv.  7^  csv. .  The  system  (2.12)  cannot  even  be  written 
in  conservation-law  form.  The  inconsistency  in  independent  variables  is  dealt  with 
directly  by  writing  a  redundant  set  of  equations  for  the  functions  s  and  cSVj . 

Transposing  x  and  y  in  (2.12a)  and  (2.12b)  results  in  the  expanded  system 


dts  +  dx 

dfCsv 


f(s)vi  +  f'{s)aSVi  +  - f"{S)(vi)0a2s 


=  0, 


dx 


f  )  CViVj  “t“  f  (s)  (c,  1  f ) C*.S‘7; 


=  0, 


dt.cs  dx 


m  Csvi  +  f(s)(Vi)0Cs 


■  a 


Vi 


f  (s'jCsvi  V  /  (^)(^)o|yCs 


=  0, 


(3.1) 


dt.cs 
dt~s  +  dVi 


dv 


f  (s^CvjVi  +  /  (S)(Uj)o|yCa 


=  0, 


f(s)vi\y  +f'(s)aSVi  +  -f”  (s)(ui)o|y5:2 


=  0; 


with  i,j  =  1,2,  where  aaVi(y,t)  =  c3Vi(y,y,t),  (f2(y ,t)  =  cs(y,y,f).  The  system  is 
symmetric  under  the  permutation  ip(pt,y,t)  =  ^>(y,x,f).  To  be  consistent  with  the 
original  system  (2.12a)-(2.12c),  the  initial  conditions  and  the  numerical  method  must 
obey  this  symmetry. 

4.  Solution.  We  numerically  solve  (3.1)  with  a  first-order  upwind  scheme  within 
the  framework  of  CLAWPACK  [21].  However,  we  do  not  take  full  advantage  of  the 
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functionality  of  CLAWPACK  clue  to  the  difficulty  in  defining  a  Riemann  solver  for 
the  MDEs.  The  solutions  are  shown  to  be  bimodal  as  in  1-D  (see  Fig.  4.1,  and  see 
[18]  for  1-D  details).  Solutions  are  also  compared  to  MCS  moment  estimates. 


Fig.  4.1.  Mean  and  standard  deviation  of  saturation  (1-D).  The  solid  curve  is  a  semi- analytical 
solution  described  in  [18];  dashed  curve  is  obtained  from  fine- grid  numerical  PDE  scheme. 


4.1.  MDEs.  Our  example  problem  is  on  a  rectangular  domain  R  =  [0,  L\]  x 
[0, L2]  =  [0^2]  x  [0,1],  x  =  (. Xi,X2 )  €  R,  with  dX2h(x)  =  0  on  aj2  =  0  and  X2  =  1 
and  constant  /i(x)  at  X\  =  0  and  X\  =  2.  The  entire  domain  R  is  initially  oil- 
saturated  (s  =  0),  with  water  (s  =  1)  pumped  in  along  x\  =  0.  We  choose  a  form 
of  the  fractional  flow  function  /(s)  that  arises  from  quadratic  relative  permeabilities: 
/(s)  =  s2/(s2  +  m(l  —  s)2).  We  set  the  viscosity  ratio  m  to  0.5  and  porosity  <j>  to 
0.2.  Statistical  parameters  for  MDEs  are  (Y)  =  0  and  a'y  =  0.25,  and  we  use  the 
exponential  covariance  function  given  by 


C(r)  =  Oy  exp 


N 

Ai 


N 

A2 


(4.1) 


where  r  =  (rq,  r2).  We  set  correlation  lengths  Ai  =  A2  =  0.2.  This  gives  5  correlation 
lengths  in  the  X2  direction  and  10  correlation  lengths  in  the  x\  direction.  Hydraulic 
head  and  velocity  moments  were  computed  using  the  semi-analytical  flow  approach  of 
Zhang  and  Winter  [30,  36] .  The  authors  derived  MDEs  using  perturbation  expansions 
in  a  manner  similar  to  ours.  Due  to  the  mathematical  complexity  caused  by  the 
presence  of  boundaries  and  nonstationarity,  these  MDEs  are  solved  numerically.  Finite 
velocity  correlation  length  is  evident  in  examples  of  velocity  covariance  shown  in  Fig. 
4.2. 

We  use  a  very  coarse  grid  (40  by  20)  due  to  the  fact  that  two-point  covariance 
functions  make  the  2-D  MDEs  essentially  4-D.  The  2-D  solutions  for  the  mean  s  and 
variance  cr2  of  saturation  are  shown  at  an  early  and  late  time  in  Fig.  4.3.  These 
indicate  a  pair  of  saturation  fronts  moving  at  distinct  speeds.  Saturation-velocity 
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Fig.  4.2.  Velocity  covariance  cvivi,  showing  finite  correlation  length.  The  solid  line  shows 
covariance  between  vi  at  the  center  of  the  domain  and  at  all  points  on  X2  =  L2/2.  The  dotted  line 
shows  covariance  between  v\  at  the  center  and  at  all  points  along  xi  =  Li/2. 


one-point  covariances  aSVi ,  i  =  1,2  are  shown  in  Fig.  4.4.  The  aSVl  surface  is  similar 
to  that  of  of,  and  is  very  close  to  the  1-D  profile. 

The  results  also  are  shown  for  s,  crSVl,  and  as  along  a  slice  at  X2  =  T2/4  on  the 
longitudinal  axis  in  Fig.  4.5.  In  each  case,  a  bimodal  solution  is  evident.  On  the  left, 
a  rarefaction  zone  trails  a  slow  front.  Physically  this  represents  a  smoothly  varying 
mixture  of  two  phases.  The  uncertainty  as  measured  by  a %  is  relatively  small  in  this 
region.  Beyond  the  slow  front,  in  the  interval  of  uncertainty,  there  is  a  nearly  constant 
ratio  of  phases.  Finally,  at  the  fast  front,  mean  saturation  and  covariances  drop  to 
zero;  beyond  the  front,  only  one  phase  (oil)  is  present. 

A  crude  test  now  can  be  performed  to  see  whether  a  covariance  relationship  holds 
similar  to  that  in  1-D.  The  flow  is  nearly  1-D;  transverse  velocities  are  much  smaller 
than  longitudinal  velocities,  and  aSV2  is  much  smaller  than  crSVl.  So  it  is  reasonable 
to  expect  that  the  products  <JgaVlVl  and  a^v  are  similar,  because  they  are  identical 
in  1-D.  However,  observe  the  dip  in  aSVl  between  two  peaks  (Fig.  4.5),  on  the  interval 
of  uncertainty.  The  saturation  variance  does  not  have  such  a  profile,  and  neither  does 
aVlVl  (see  Appendix  A).  Thus  there  is  a  difference  in  shape  between  the  two  products 
<j'j<JVlVi  and  <TgVi,  so  that  the  1-D  identity  (T^aVlVl  =  cr1SVi  cannot  hold. 

The  2-D  solutions  are  comparable  to  solutions  in  the  limit  of  1-D  flow.  This 
limiting  solution  is  obtained  by  letting  the  longitudinal  velocity  components  V\ ,  (t>i)o 
be  constant  and  setting  transverse  velocities  to  zero.  The  limiting  solution  matches 
a  1-D  solution  computed  on  the  coarse  grid.  The  result  along  a  longitudinal  slice 
is  shown  in  Fig.  4.6.  (The  1-D  result  is  diffused  compared  to  Fig.  4.1  because  of 
the  coarse  grid.)  The  nearly  uniform  2-D  flow  scenario  chosen  gives  results  nearly 
identical  to  those  in  1-D  in  spite  of  the  localized  velocity  correlation  structure  in  2-D. 

The  addition  of  small  linear  diffusion  terms,  which  are  employed  in  the  numerical 
scheme  to  stabilize  the  solution,  does  not  eliminate  bimodality.  The  fact  that  we  solve 
on  a  coarse  grid,  introducing  large  numerical  diffusion,  makes  it  even  more  striking 
that  bimodal  solutions  are  clearly  evident  in  the  solution.  Grid  refinement  and  more 
accurate  numerical  methods  can  be  expected  to  emphasize  the  wave  separation  even 


more. 
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(a)  Mean  saturation 
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Fig.  4.3.  Saturation  mean  and  variance  contours  at  (a,b)  early  and  (c,d)  late  times.  Lighter 
shades  indicate  a  higher  mean  water  saturation  or  higher  variance. 


The  separation  between  wave  fronts  increases  with  time.  In  many  specific  appli¬ 
cations  the  relevant  time  scale  may  be  such  that  a  separation  is  not  evident;  however, 
even  when  diffusion  dominates  initially,  bimodality  will  appear  over  sufficiently  long 
times.  This  property  of  the  solution  must  be  accounted  for  in  any  use  of  MDEs  for 
advection  and  diffusion,  for  both  single-phase  and  multiphase  applications. 

4.2.  Comparison  to  MCS.  Monte  Carlo  simulations  are  frequently  used  in  ap¬ 
plications  to  subsurface  flow  and  transport  and  for  reservoir  modeling.  The  primary 
drawback  is  that  MCS  is  computationally  costly.  Separate  and  very  important  issues 
include  the  problem  of  fair  sampling  of  the  true  sample  space  and  the  standard  prob¬ 
lem  of  autocorrelation  in  random  number  generators.  However,  it  is  a  well-established 
method  and  serves  as  a  standard  for  comparison. 

For  a  bound  on  the  error  |s(x,  t)  —  ( s )  (x,  f)|  of  just  10-2,  we  estimate  that  about 
4000  simulations  are  needed.  This  is  a  significant  burden  on  time  and  space  resources, 
so  we  present  results  of  only  350  simulations.  This  gives  an  estimated  absolute  error 
bound  of  about  0.032,  which  is  sufficient  for  a  rough  comparison  (recall  that  mean 
saturation  values  are  restricted  to  [0, 1]). 

Local  correlation  structure  is  again  evident  in  MCS  velocity  moments  in  Fig.  4.7 
(compare  to  Fig.  4.2).  Sample  mean  and  an  envelope  of  one  standard  deviation  along 
the  midline  aq  =  L2/2  are  shown  in  Fig.  4.8.  Contours  of  sample  moments  are  shown 
in  Fig.  4.9.  The  roughness  of  the  contours  highlights  the  fact  that  a  greater  number 
of  simulations  is  required  to  obtain  high  accuracy.  The  mean  front  is  near  aq  =  1.3, 
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Fig.  4.6.  Comparison  of  2-D  and  1-D  flow. 


X./L. 

I  l 


Fig.  4.7.  Sample  velocity  covariance  CV1V1,  showing  finite  correlation  length.  The  solid  line 
shows  covariance  between  tii  at  the  center  of  the  domain  and  at  all  points  along  x,2  =  L2/2.  The 
dotted  line  shows  covariance  between  vi  at  the  center  and  at  all  points  along  x\  =  Li/2. 


where  the  variance  is  greatest. 

To  make  the  comparison  between  MCS  and  MDE  as  consistent  as  possible,  the 
MDE  solution  is  computed  here  using  the  sample  velocity  moments.  This  comparison 
isolates  the  difference  between  the  two  methods  at  the  level  of  computing  saturation 
moments  (kindly  suggested  by  Neuman  and  Guadagnini  [24]).  Contours  of  MDE 
solutions  using  MCS  velocity  moments  are  shown  in  Fig.  4.9,  where  they  are  compared 
to  MCS  saturation  moments.  A  clearer  comparison  is  seen  in  profiles  of  saturation 
moments  in  Fig.  4.10,  which  shows  solutions  at  a  late  time  along  X2  =  T2/2.  Of 
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Fig.  4.8.  Saturation  sample  mean  profile,  with  an  envelope  of  one  standard  deviation. 


Fig.  4.9.  MCS  saturation  (a)  mean  and  (b)  variance  and  MDE  (c)  mean  and  (d)  variance 
using  MCS  velocity  moments. 


course,  MCS  and  our  MDE  cannot  match;  there  is  no  physical  reason  to  expect  that  a 
bimodal  solution  can  be  obtained  from  simulations  in  our  model.  However,  the  MDE 
mean  is  centered  near  the  simulation  mean,  and  variances  computed  by  both  methods 
are  similar.  The  MDE  solution  is  ahead  of  the  MCS  solution,  which  is  likely  due  to 
the  effects  of  numerical  and  artificial  diffusion.  The  first  of  these  effects  is  much  larger 
for  MDE  due  to  the  coarse  grid  on  which  it  is  computed. 

The  comparison  shows,  not  surprisingly,  that  the  MDE  moments  are  not  an  ac¬ 
curate  representation  of  the  details  of  MCS  moments.  They  provide  a  useful  first 
approximation,  particularly  if  a  rough  estimate  of  the  location  and  magnitude  of  un¬ 
certainty  is  sufficient  in  an  application.  The  comparison  is  favorable  in  the  rarefaction 
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Fig.  4.10.  Comparison  of  MDE  and  MCS  saturation  moments. 


zone,  but  we  are  more  concerned  with  behavior  near  the  (true)  mean  saturation  front. 
In  1-D,  [18]  explains  the  observed  bimodality  in  terms  of  the  eigenvalues  and  eigen¬ 
vectors  of  the  hyperbolic  system  of  MDEs.  This  finite-dimensionality  results  from  the 
truncation  of  the  expansions  in  (2.3)  and  (2.4).  We  speculate  that  if  more  terms  were 
retained,  the  number  of  modes  would  increase,  but  their  magnitude  would  decrease, 
forming  a  closer  representation  of  MCS  moments. 

For  integrated  or  averaged  quantities,  the  MDE  estimates  are  better.  For  example, 
the  solution  does  provide  a  good  estimate  of  the  total  production,  which  is  a  key 
measure  in  the  oil  industry.  The  total  oil  produced  is  computed  as  the  mass  quantity 
that  has  exited  the  right  boundary  at  time  t.  The  volume  of  oil  this  represents  is  the 
same  as  the  volume  of  water  injected  in  this  case,  which  is  given  by 

P{t)  =  JJ  <f>s(x,t)dA,  (4.2) 

(assuming  unit  domain  length  in  the  vertical  and  unit  density  of  oil)  in  each  simulation. 
Sample  mean  and  variance  of  P(t)  are  computed  and  compared  to  an  estimate  of  mean 
production  from  MDE.  The  latter  is  obtained  by  replacing  s  in  the  integral  with  s. 
The  results  are  shown  in  Fig.  4.11  (here  only  200  simulations  were  used). 

The  curves  are  fairly  closely  matched;  the  MDE  curve  lies  within  one  sample 
standard  deviation  ds  of  the  MCS  curve,  except  near  t  —  0.  The  MDE  curve  has 
two  sharp  changes  in  slope,  which  are  due  to  the  passing  of  the  two  saturation  fronts 
through  the  right  boundary.  In  each  of  the  first  two  segments,  the  curve  is  nearly 
linear,  which  it  should  be  for  an  intransient  velocity  field.  Until  water  reaches  the 
boundary  (breakthrough),  a  constant  flux  of  oil  leaves  the  domain.  The  flux  is  again 
nearly  constant  between  the  saturation  fronts,  providing  the  second  linear  segment. 
Finally,  after  the  second  front  passes  the  boundary,  production  decays  smoothly  as 
the  rarefaction  part  of  the  mean  saturation  crosses. 
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Time 


Fig.  4.11.  Total  Production.  Solid  curve  is  obtained  by  MDE,  dashed  curve  is  MCS  mean.  One 
standard  deviation  is  shown  by  dotted  curves  (200  samples). 


MCS  provides  solutions  with  high  resolution;  this  is  also  a  requirement,  as  the 
velocity  field  in  each  simulation  may  change  significantly  over  short  spatial  scales. 
Averaging  provides  relatively  smooth  velocity  and  saturation  moments.  The  MDE 
approach  takes  advantage  of  the  smoothness  of  these  moments  before  they  are  com¬ 
puted.  Thus,  we  do  not  need  a  fine  grid  to  capture  detail.  Grid  refinement  will 
certainly  improve  the  solution  to  MDE,  but  refinement  is  not  necessary  to  capture 
small-scale  features. 

We  computed  rough  estimates  of  computational  costs  and  found  that  on  an  SGI 
300MHz  MIPS12K  processor  with  256M  memory,  a  set  of  200  simulations  on  200 
by  100  grid  nodes  takes  about  8  hours.  To  achieve  an  error  tolerance  of  ICC2  (to 
within  97.5%  confidence),  we  estimate  that  4000  simulations  are  necessary,  which 
would  take  about  160  hours  on  the  same  computer.  This  is  much  longer  than  the 
5  hours  it  takes  to  obtain  the  solutions  presented  here  using  the  MDE  code.  These 
are  crude  estimates;  a  careful  comparison  of  costs  for  the  two  methods  must  allow  for 
some  attention  to  improving  efficiency  in  both  codes  and  should  account  for  several 
differences  between  the  methods,  including  the  significantly  higher  storage  cost  for 
MDE.  Additional  details  on  the  MCS  approach  we  used  can  be  found  in  Appendix  B 
and  in  [16]. 

5.  Conclusions.  Bimodal  moments  are  obtained  for  saturation  in  2-D  as  in 
earlier  results  for  1-D,  in  spite  of  the  finite  length  scale  for  velocity  correlations  in 
2-D.  The  fact  that  velocity  correlations  lead  to  a  macrodispersion  term  in  lower-order 
approximations  to  moments  would  suggest  a  partial  transition  from  hyperbolic  to 
parabolic  character  as  we  move  from  1-D  to  2-D.  However,  in  the  example  herein 
we  limit  the  domain  to  ten  correlation  lengths  (far  from  the  ergodic  regime)  so  that 
macrodispersion  effects  may  be  very  limited. 

The  two-phase  Eulerian  MDEs  presented  here  complement  the  traditional  MCS 
approach.  MDEs  provide  a  first  approximation  to  moments  relatively  quickly.  Rather 
than  compute  several  hundred  sample  saturation  fields  on  fine  grids  and  post-process 
to  obtain  moments,  a  single  solution  of  MDEs  is  computed  on  a  coarse  grid.  In  spite  of 
a  non-physical  bimodality  in  MDE  moments,  we  obtain  a  good  match  between  MDE 
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and  MCS  in  the  total  oil  production  curve.  For  high  accuracy,  MCS  is  appropriate 
and,  in  real  applications,  it  may  be  of  most  benefit  to  combine  the  two  methods.  It 
may  also  be  possible  to  improve  the  accuracy  of  MDEs  by  retaining  more  terms  in 
the  perturbation  expansions. 

The  Eulerian  MDE  approach  applies  to  any  probability  distribution  of  geologic 
properties  with  any  correlation  function.  It  depends  solely  on  moments  rather  than 
a  particular  probability  distribution  and  does  not  require  stationarity,  thus  removing 
restrictions  that  are  common  to  stochastic  subsurface  theories. 

6.  Acknowledgments.  We  are  grateful  to  Joseph  Oliveira  for  many  conversa¬ 
tions  on  this  topic. 

Appendix  A.  Velocity  moments  from  MDE. 

Additional  velocity  moment  profiles  are  shown  in  Fig.  A.l  and  Fig.  A. 2.  These 
were  obtained  using  the  MDE  flux  code  of  Zhang  and  Winter  [36].  The  dominant 
component  is  longitudinal.  Transverse  flow  is  inward  near  the  inflow  (left)  boundary, 
and  outward  near  the  outflow  (right)  boundary.  Longitudinal  flow  is  significantly 
greater  along  the  sides,  where  transverse  flux  is  zero. 


(a)  (b) 
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Fig.  A.l.  Mean  velocity  estimates  ( t?o)i  and  v\  at  (a)  x 2  —  L2/2  and  (b)  x\  —  L\/2;  velocity 
variances  crViVi  at  (c)  X2  =  £2/2  and  (d)  x\  =  L\/2. 
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x  10'3 


Fig.  A. 2.  Mean  transverse  velocity. 


Appendix  B.  MCS  and  convergence. 

Monte  Carlo  simulations  were  performed  by  numerically  solving  the  total  velocity 
and  saturation  equations  (1.1)  and  (1.2)  over  replicates  of  log  hydraulic  conductiv¬ 
ity  with  correlation  function  (4.1).  Replicates  were  generated  with  a  Fast  Fourier 
Transform  algorithm  developed  by  Gutjahr  and  colleagues  [14].  Velocity  equations 
were  solved  with  a  conjugate  gradient  algorithm.  Saturation  profiles  were  obtained 
using  a  first-order  upwind  approach  just  as  for  MDEs,  also  within  the  framework  of 
CLAWPACK. 

We  performed  simple  tests  for  convergence  of  the  mean  and  variance  of  saturation. 
A  tolerance  of  approximately  0.032  to  within  97.5%  confidence  is  estimated  for  results 
of  the  350  simulations  presented  here.  Cumulative  average  saturation  from  the  first 
150  simulations  at  representative  grid  points  is  shown  in  Fig.  B.l.  Although  it  is  well- 
known  that  convergence  of  such  estimates  may  appear  imminent  when  in  fact  they 
are  not,  this  is  still  a  useful  tool.  We  may  be  reasonably  confident  in  the  simulation 
moments  for  comparison  to  MDE;  however,  it  is  clear  from  the  figures  in  §4.2  that 
a  larger  number  of  simulations  would  be  needed  to  obtain  the  smooth  profiles  we 
expect. 


(a) 


(b) 


Fig.  B.l.  Convergence  of  sample  mean  saturation  at  (a)  X2  =  0.4  and  (b)  X2  =  0.4,-  in  both, 
xi  =  1.2  (solid),  xi  =  0.95  (dashed),  and  x\  =  0.7  (dotted). 
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